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Abstract 

This paper introduces the FermiFab toolbox for many-particle quantum systems. It 
is mainly concerned with the representation of (symbolic) fermionic wavefunctions and 
the calculation of corresponding reduced density matrices (RDMs). The toolbox trans- 
parently handles the inherent antisymmetrization of wavefunctions and incorporates the 
creation/annihilation formalism. Thus, it aims at providing a solid base for a broad 
audience to use fermionic wavefunctions with the same ease as matrices in Matlab, say. 
Leveraging symbolic computation, the toolbox can greatly simply tedious pen-and-paper 
calculations for concrete quantum mechanical systems, and serves as "sandbox" for the- 
oretical hypothesis testing. FermiFab (including full source code) is freely available as a 
plugin for both Matlab and Mathematica. 
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1. Introduction 

The ground state energy of fermionic many-particle quantum systems can be re- 
expressed as a linear functional of (one- or two-body) reduced density matrices (RDMs) . 
This notion traces back to the origins of quantum mechanics [1, 2] around 1930. Since 
1964, the one-body RDM has been greatly popularized by density functional theory [3, 4], 
which is typically the most viable approximation for handling large particle numbers. 
The tantalizing possibility of employing RDMs (instead of many-particle wavefunctions) 
for exact groundstate energy computations is counterbalanced by the iV-representability 
problem, i.e., the search for necessary and sufficient conditions a two-body density must 
obey to represent an A'^-electron wavefunction [5, 6, 7]. Modern applications use varia- 
tional principles and semidefinite programming to impose positivity constraints on the 
two-body RDM [8]. In any case, it is desirable to render the powerful RDM framework 
accessible to a broader audience, integrating it into the symbolic language of modern 
computer algebra systems like Mathematica, or numeric software like Matlab. 

The FermiFab toolbox (available for download at [9]) is precisely designed for that 
purpose. A short "usage manual" and a brief tour of the essential features is provided 
in the following subsections. Note that the underlying one-particle orbitals (see below) 
are always assumed to be orthonormalized. In addition, the toolbox adheres to the 
trace-normalization convention tr^p-^7|^)^^| = (^) for the p-body RDM 7|^)(^| of a 
normalized A^-body wavefunction i/j. Here, A^H denotes the p-particle Fock-space (see 
following subsection). 

1.1. Fermi states 

Fundamental building blocks of multi-fermion quantum systems are Slater determi- 
nants (figure 1). These can be thought of as a collection of "orbitals" (or slots), some 
of which are occupied by a fermionic particle (e.g., an electron). In mathematical terms, 
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Figure 1: Schematic illustration of a Slater determinant: (filled) circles correspond to (occupied) orbitals. 

the available number of orbitals 'orbs' is the dimension of the underlying one-particle 
Hilbert space H and the number of occupied orbitals the particle number N. Thus there 
are altogether {°'^^^) Slater determinants. Their complex span defines the A^-particle 
Fock-space A^H. The A^-particle Fermi states are precisely the elements of A^H. 

1.2. Getting started with FermiFab 

For concreteness, the following examples are issued in the Matlab programming lan- 
guage. (The Mathematica version of FermiFab provides the same features; section 3 
contains a demonstration.) Commands typed by the user are preceded by », and the 
subsequent lines show the corresponding output. In standard Matlab syntax, zeros (n, 1) 
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below constructs a column vector of length n, and nchoosek computes binomial coeffi- 
cients. We first show how to represent an A'' = 4 particle state tp with, e.g., 6 available 
orbitals in total: 

» orbs =6; N = 4; 

» X = zeros (nchoosek (orbs, N) ,1) ; x(l)=l/sqrt (2) ; x(2)=li/sqrt(2) ; 
» psi = f ermistate(orbs,N,x) 

psi = 

Fermi State (orbs ==6, N == 4) 
(0.70711) I 1234> + (0+0.707111) I 1235> 

Needless to say, the fermistate command is specific to the FermiFab toolbox. The 
vector X contains the Slater determinant coefficients of in lexicographical order. Let's 
assign more meaningful names to the orbitals of ip: 

» psi = set(psi, 'orbnames' ,{'ls' 'ls~' '2s' '2s~' '2p' '2p~'}) 
psi = 

Fermi State (orbs ==6, N == 4) 

(0.70711) Us Is" 2s 2s~> + (0+0 . 707111) I Is ls~ 2s 2p> 

From a physics viewpoint, these orbitals could form electronic subshells in atoms. The 
rank-one projector \tl)) {ip\ or "density matrix" of ip can be calculated intuitively by 

» psi*psi' 
ans = 

Fermi Operator weclge"4 H -> weclge"4 H (orbs == 6) 

Matrix representation w.r.t. ordered Slater basis 

(Us Is" 2s 2s~>, Us ls~ 2s 2p>, ... 1 2s 23" 2p 2p~>) -> 
(Us ls~ 2s 2s''>, Us ls~ 2s 2p>, ... |2s 2s" 2p 2p~>) : 

Columns 1 through. 4 

0.5000 - 0. 50001 



Note that the result is now a f ermiop operator acting on A H. 

1.3. Reduced density matrices 

The core feature of the toolbox is the efficient calculation of RDMs. For example, the 
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2-body RDM 




can be obtained by 



» rdm(psi,2) 



ans = 
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Fermi Operator wedge~2 H -> wedge"2 H (orbs == 6) 



Matrix representation w.r.t. ordered Slater basis 

(Us ls~>, lis 2s>, ... |2p 2p~>) -> (lis ls''>, lis 2s>, 



|2p 2p->): 



Columns 1 through 4 



1.0000 









1 . 0000 












0.5000 
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RDMs are reviewed in more detail in section 2.3. 

Tensor products of operators 

Given a linear operator A : H — "H, a straightforward derivation based on the 
antisymmetrized structure of A^H shows that 



for all 1 < ii < • • • < z^v < dimH and 1 < ii < • • • < < dimH. That is, we 
obtain a matrix representation of ■ ■ ■ ^ A acting on A^H. The tensor _op command 
implements precisely this operation. The following code lines are taken from the "natural 
orbitals" example in test/norbs .m: 

» orbs =6; N = 4; 

>> psi = f ermistate (orbs ,N, crand(nchoosek(orbs ,N) , 1) ) ; 
» [U,D] = elg(rdm(psi,l)) ; 

crand generates pseudorandom complex numbers (similar to rand), and eig computes 
eigenvalues and -vectors. Thus, the eigenvectors of the 1-body RDM of ip are stored in 
U. Performing a corresponding base change on A^T-l using these eigenvectors should 
result in a diagonal 1-body RDM [5] : 

» psi = tensor_op(U,N) '*psl; 
» G = get (rdm(psl, 1) . 'data') ; 
>> err = norm(G-dlag(dlag(G) ) ) 



1.6512e-015 

In many physical applications, one can take advantage of unitary base changes on H such 
that subsequent computations are simplified, e.g., by choosing single-particle eigenstates 
of the Lz angular momentum operator. The above code shows how to implement the 
according base change on A^T-L. 



{ji, . . ■ ,3n\{A® ■ ■ ■ ® A)\ii, . . . ,iN) 



det Ofc M|«f)fc,^ 



err = 
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1.5. State configurations 

For performance and memory efficiency reasons, FermiFab has built-in "configura- 
tions", i.e., we can subdivide the available orbitals into several groups, each of which 
contains a fixed number of particles. (Physically speaking, the groups could be inter- 
preted as atomic subshells Is, 2s, 2p, 3s, for example.) Let's say our system involves a 
total of 3 particles in 9 orbitals, with exactly 2 particles in the first 5 orbitals and 1 
particle in the remaining 4 orbitals. Then a f ermistate refiecting this configuration is 
specified by 

» orbs = [5,4] ; N = [2,1] ; 
» psi = f ermistate (orbs, N) 

psi = 

Fermi State (orbs ==9, N == 3) 
|126> 

Note that |126) is the lexicographically first base vector respecting the configuration 
constraints, and that V' requires only (j) • (^) = 40 rather than (g) = 84 complex entries: 

» length. (psi) 
ans = 
40 

The rdm command works transparently for any configuration, so ij) behaves like a standard 
9-orbital, 3-particle state. 

What happens if we add two f ermistates with different but compatible configura- 
tions (i.e., the total number of orbitals and particles is the same)? 

» orbs = [2,7] ; N = [1,2] ; 
» phi = f ermistate (orbs, N) 

phi = 

Fermi State (orbs ==9, N == 3) 
|134> 

» length (phi) 
ans = 
42 

» chi = psi+phi 
chi = 

Fermi State (orbs ==9, N == 3) 
|126> + |134> 

as expected so how is this accomplished? FermiFab has detected that it needs to 
combine the two configurations, resulting in the full-fiedged 9-orbital, 3-particle state. 
This fact can be checked by 
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» length ( Chi) 
ans = 
84 

1.6. Symbolic computations 

The Mathematica version of FenniFab is - quite naturally - inherently based on sym- 
bolic language. Considering Matlab, the (optionally available) Symbolic Math Toolbox 
integrates seamlessly into FermiFab, too. Taking advantage of symbolic computations is 
thus easily accomplished. That is, in the above examples, we may as well insert symbolic 
variables: 

» syms a b c 

» y = symCzeros (1 ,nchoosek(orbs ,H) ) ) ; 
» yd) = a; y(3) = li*b-2; y(4) = 1/c; 
» psi = setCpsi, 'data' ,y) 

psi = 

Fermi State (orbs ==6, N == 4) 

(a) I Is Is" 2s 2s~> + (b"2*i)|ls is" 2s 2p"> + (1/c) I Is is" 2s" 2p> 
» rcl]n(psi,2) 
ans = 

Fermi Operator wedge"2 H -> wedge"2 H (orbs == 6) 

Matrix representation w.r.t. ordered Slater basis 
(|12>, |13>, ... |56» -> (|12>, |13>, ... |56»: 

[ (c*b*2*conj (b)*2 + a*c*conj (a))/c + l/(c*conj (c)) , 



2. Implementation Details 

The algorithmic implementation is based on the canonical mapping from Slater de- 
terminants to bitfields. That is, each Slater determinant corresponds to an unsigned 
integer s, where the ith bit is set to 1 precisely when the ith orbital is occupied. To 
remain unambiguous in terms of bitlength, the first orbital is stored in the LSB (least 
significant bit). Now, our task consists of re-expressing the creation/annihilation and 
RDM formalism in terms of bit operations. Note that, for example, testing whether all 
occupied orbitals in si are also occupied in s^ amounts to the pretty simple line of code 
si A S2 = si, where we have used the bitwise AND operator A. The following table 
summarizes all required bit operations: 
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bitwise AND: x ^y 

bitwise OR: xW y 

bitwise XOR: x®y 

bit shift left: x ^ n 

bit shift right: x ^ n 

bit count: #(2;) 



For example, IOOIIOI2 > 3 = IOOI2 and #(18) ^ #(100102) = 2. Note that bit 
operations are typically very "cheap" on CPUs. (In particular, refer to the SSE4 [10] 
POPCNT "population count" instruction for bit counting.) Diving a little bit further 
down into CPU intrinsics, we will make use of two's-complement arithmetic for negating 
numbers [11], e.g., 

...OOIOIIIOO2 
-x = ... IIOIOOIOO2. 

Interestingly, precisely all bits flip which are more significant than the least significant 
1-bit (marked red). Thus, we can use this property to extract the last 1-bit from a bitfield 
X ^ Q simply by 

LastBit(a::) := x A {~x). 

(An less universal alternative is the BSF "bit scan forward" instruction [12], which returns 
the index of the least significant 1-bit.) 

2.1. Enumerating Slater determinants 

The basic task we set out to accomplish in this subsection is lexicographically enu- 
merating all Slater determinants of a fixed particle number N and number of orbitals 
'orbs'. This amounts to computing the lexicographically next bit permutation (denoted 
by 'NextFermi'). For example, 

s = OIIIIOOO2 
NextFermi(s) = 100001 II2. 

Closer inspection reveals the general rule that the leading 1-bit (marked red) in the 
least significant block of Is gets shifted to the left by one position, and the remaining 
1-bits are shifted to the end. Algorithm 1 is adopted from [13] and performs exactly 
this computation. In line 1, s V (s — 1) sets the trailing zeros in s to 1, so for example, 
s = OIIIIOOO2 s V (s - 1) = OIIIIIII2 and t = IOOOOOOO2. The second term in line 2 
performs the shifting of the remaining 1-bits to the end. 



Algorithm 1 NextFermi 
Input: s: bitfield 

1: t ^ (S V (S - 1)) + 1 

2: return t V (((LastBit(t) - l)/LastBit(s)) > 1) 



As an extension of Algorithm 1, we want to take into account "configurations", i.e., 
a subdivision of the available orbitals into several groups, each of which contains a fixed 
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number of particles. For example, we compartmentalize a total of 11 orbitals such that 
exactly 4 particles are in the first 6 orbitals and 2 in the remaining 5 orbitals, written as 
(orbsi, orbs2) = (6,5) and (iVi, A^2) = (4,2). Then a sequence of patterns - respecting 
the configuration restrictions - would be 

0|01010|1101102, 
0|01010|11100l2, 

0|01010|1110102, (2) 

0|01010|1111002, 

0|01100|00111l2, 

where we have highlighted the currently changing 1-bits by red colors. 

More formally, given (orbsi, . . . ,orbsfc), the compartmentalization may be written as 
Vj := span{|i) : bj-i < i < bj} C H with bj := J2e=i orbsf. In other words, H — 0^ Vj. 
In the example above, Vi = span{|l) , . . . , |6)} and V2 = span{|7) , . . . , |11)}. Now, 
mathematically speaking, a configuration of an A^-particle state is a subspace of A^H of 
the following form: 

CNu...,N, :=span{|zi,...,ZAr) : ■ & V,} = N^} (3) 

where (iVi, . . . ,Nk) is a partition of N (i.e. < Nj < orbsj, Nj = N). A quantum 
chemist could interpret the Vj as atomic subshells Is, 2s, 2p, 3s, . . . and the Nj as occu- 
pation numbers. An interesting consequence of definition (3) is the recovery of a tensor 
product structure, namely 

k 

This follows from the observation that a configuration is constructed by the lexicograph- 
ical enumeration of Slater determinants within orbital groups, as illustrated in (2). 

Algorithm 2 implements precisely this enumeration. In accordance with the lexico- 
graphical scheme, it first iterates through all Slater determinants within the least signif- 
icant orbital group (line 3), then resets this group (first term in line 12) and recursively 
computes the next bit pattern for the remaining groups (line 8). The mask in line 1 is 
required for testing whether the last bit permutation within the least significant group 
has been reached (line 2). In the example above, we would have mask = 0|llllll2. 

2.2. Creation/annihilation operators 

The creation/annihilation operator formalism is an essential ingredient of many- 
particle quantum mechanics and quantum field theory [14]. For a very brief sketch, 
let if S A^Ji be a p-particle wavefunction with 1 < p < N . Then, the linear annihila- 
tion operator acting on A^Ti, removes or "annihilates" the state ip from A^T-L. More 
precisely, a^p is uniquely determined by its antilinearity in Lp, 

Vce C,(^i,v?2 e Am 

together with the decomposition for Slater determinants. 



a\ii,i2,...ip) ■= a|ip) ■ ■ ■a\i2)0'\ii) V 1 < ii < • • • < ip < dimH, 
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Algorithm 2 NcxtFermiConfig 



Input: s: bitfield, orbs: int array 
1. mask (1 < orbs[0]) - 1 
2: if (((s V (s — 1)) A mask) 7^ mask) then 
3: return NextFermi(s) 
4: else 

5: if orbs. length — 1 then 
61 return -1 
7: end if 

8, t <^=: NextFermiConfig(s orbs [0], orbs [1, . . . , end]) 

9: if t = -1 then 
10: return -1 
11: end if 

12: return (mask/LastBit(s)) V (i ^ orbs[0]) 
13: end if 



as well as the definition 

for all 1 < < • • • < < dim'H. The sign factor can be interpreted as number of 
orbital "flips" illustrated in figure 2. 



(-1) (-1) (-1) 



6 7 
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Figure 2: Annihilation of a single orbital (red). Figuratively, the red orbital moves to the front before 
being removed, such that each flip (curved arrows) with an occupied preceding orbital contributes a sign 
factor of (—1). In terms of quantum mechanics, ajg^ |24568) = — |2458). 

So far we have considered annihilation operators only. The creation operator is 
by definition the adjoint (conjugate transpose) of a^p^ as the notation already suggests. 
It can be shown that the following relations hold, where the anticommutator bracket is 
defined by {A, B] := A B — B A and (p,x ^ A^H are arbitrary wavefunctions: 

{a^,a^}^0, {al,al}^0, {a^, aj^} = ((^ | x) ■ 

In the remainder of this subsection, we want to detail an efficient algorithmic im- 
plementation of the annihilation operation, w.l.o.g. for Slater determinants only. More 
precisely, let \s) e be a fixed Slater determinant, then our task is the calculation 

of a|j) \s) for arbitrary Slater determinants \t) G A^H and 1 < p < N. The result will be 
nonzero only if all occupied orbitals in t are also occupied in s, which can be tested by 
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t A s = t a.s already mentioned in the beginning. Given this holds true, the bit pattern 
describing the Slater determinant a\t^ \s) is simply s — t, so what essentially remains is 
the calculation of the sign factor. 

For that purpose, we define the annihilation sign mask of s such that each bit stores 
the integer parity of the number of less or equal significant 1-bits in s. That is, if s has 
binary representation 

s =...020100, OiS{0, 1}, then 

i 

AnnihilSignMask(s) := . . . 626160 with 6^ = aj mod 2. 

i=o 

76543210 76.543210 

For example, s =001011002 results in AnnihilSignMask(s) = . . . IIIOOIOO2 where blue 
overhead numbers label bit positions. Algorithm 3 implements this calculation. It has a 
running time of C'(#(s)) since the last statement (line 5) in the while loop removes the 
least significant 1-bit from s. In line 4, the (B{—t) operation flips all bits which are less 
or equal significant than the current least significant 1-bit. 



Algorithm 3 AnnihilSignMask 
Input: s: bitfield 

1: bitfield m <J= 

2: while s 7^ do 

3: t <= LastBit(s) 
4: m <^= TO ® {—t) 

5: S ^ S — t 

6: end while 
7: return to 



Finally, we define the reverse permutation sign CTrcvperm for a-U ^ G N>i by the 
sign of the permutation i 1-^ n — i + 1 (i = l,...,n). A moment's thought reveals that 

fT (n) - f-1l5("-i)" 

'-'revperm V / — V / 

Altogether, our devised algorithm is illustrated in figure 3. More formally, we obtain 
'^{t) \s) — C ' \s — t) with the sign factor C equal to 

C = ^rcvpenn (#(t)) • (-l)*<'''--^*\ (5) 

where we have set Omask '■— AnnihilSignMask(s) <C 1. Equation (5) will be the basic 
building block for calculating reduced density matrices in Algorithm 4 below, as described 
in the next subsection. 

2.3. Reduced density matrices 

In this subsection we briefly recall the relevant abstract formalism, and then describe 
the algorithmic implementation in the FermiFab toolbox. Let 1 < Pfe < (fc = 1,2) 
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Figure 3: The cumulative sign factors incurred during rearrangement of the to-be annihilated (red) 
orbitals to the front. In terms of quantum mechanics, the corresponding operation reads a|45g) |24568) = 
1 26). The contribution from all flips (curved arrows) during each step can be obtained from the marked 
bit in ctmask- Note that this mask needs to be calculated only once. The permutation sign for sorting the 
three red orbitals in the last step equals cTrcvperm (3) = —1, so the overall sign factor is 1. Algorithmically, 
the whole schematic is implemented by equation (5). 

and denote orthonormal basis sets of AP''H by (ipki)i- For wavefunctions tpk G A^'='H 
(fc = 1,2), define the reduced density matrix 7|^i)(^2| : A^^H — > A^^T-L by 



(6) 
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where we have employed the creation/annihilation operators defined in the last sub- 
section. The significance of this definition can be seen as follows. Any linear map 
b : A^^H K^^T-L with matrix representation (6^) may be "lifted" to an operator 
B : A^i-H ^ by 

B --^YlKal^^a^,^. (7) 

(A prominent example is the Coulomb operator {pi = P2 = 2), which describes the 
pairwise interaction between charged particles.) Now, the B expectation value with 
respect to (^'2! equals 

In other words, this equation switches from A^'^'H to A^'^H {k = 1,2). For many appli- 
cations, this is the only possibility to avoid the "curse of dimensionality" induced by the 
Ni, Af2-particle systems. In terms of FermiFab, (7) is implemented by the p2N command. 

In the rest of this subsection, we focus on the calculation of 7|^i)(^2| in Algorithm 4. 
Due to linearity, it suffices to restrict ourselves to Slater determinants. That is, ipi and ip2 
are (w.l.o.g.) replaced by Slater determinants si and S2, respectively, and it is assumed 
that the (fki) are Slater determinants, too. So the last term in (6) can be concisely 
written as {at^S2 \ a-tiSi) with Slater determinants G A^'^H {k — 1,2). Note that the 
particle number conservation law imposes Ni—pi — N2 — P2 , otherwise all terms will be 
zero; so we calculate p2 from given A^i, N2 and pi. 



force choice 




Figure 4: Alignment of two Slater determinants for the annihilation operation. "Force" labels the orbitals 
which are either occupied in si or S2, but not in both, whereas "choice" labels all orbitals occupied in 
si as well as S2- 

The basic algorithmic idea is exemplified in figure 4. Namely, we subsume all orbitals 
occupied either in si or S2, but not in both, as "force" group, and all orbitals occupied 
in both si and S2 as "choice" group. The corresponding bit patterns /mask and Cmask are 
computed in lines 1 and 13 of Algorithm 4 by a single bit operation. Since {at2S2 \ o-tiSi) 
is nonzero only if a^jSi = ±0(332, all occupied "force" orbitals have to be annihilated by 

and at^ , respectively. On the other hand, each "choice" orbital annihilated by at^ 
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must also be annihilated by at^ and vice versa, but there's a freedom in exactly which 
of these orbitals to select, hence the "choice" designator. In our example, the only force 
orbital occupied in si is 7, so ti must contain 7 but may "choose" between 5, 6 and 9. 
If Pi = 3, we obtain ti equal to one of |567), |579) or |679). The respective t2 states 
are then |1245 6), |1245 9) and 1 1 2 4 6 9) . After the obligatory annihilation sign factor 
calculations (5), the final result (for pi — 3) is 

715679X1245691 = 1567) (124561 -1579) (124591 -1679) (12469|. 

Algorithm 4 implements equation (5) in line 9 and the first term of line 19. SforcG,fc 
stores the orbitals which must be annihilated in Sk {k — 1,2), and the number of to-be 
annihilated "choice" orbitals in si is computed in line 3. The while loop accumulates 
the return value list r containing the ket-bra's as in the above example. In line 18, 
the algorithm uses the 'BitDistribute' command, which basically just shifts bits to the 
positions designated by the 1-bits in Cmask- 

Algorithm 4 SlaterRDM 
Input: si, bitfield, pi: int 

1: /mask Si ® S2 // "forcc" mask 

2: Storcc,fc (/mask A Sfc) (fc = 1, 2) 
3: Tlchoice,! Pi — #(sforco,l) 

4: if rtchoico,! < then 
5: return 
6: end if 

7: P2-^#(S2)-#(S1)+Pl 

8: Oniask.fc AnnihilSignMask(sfc ) ^ 1, A; = 1, 2 

nLl f^rcvpcrm (Pk) ' (- l)#(''»-k,. Asf„ec..) // ^ign factor 

10: if nchoicc.i = then 

11: return C ' kforco.l) (Sforcc,2| 

12: end if 

Cmask ^ Si A S2 // "choice" mask 

14: /^choice ^(Cmask) 

15: r ^ {} 

16: i (K 71choicC,l) - 1 

17: while (t ^ fcchoicc) = do // iterate Fermi map of 'choice' orbitals 
18: Schoico BitDistribute (t, Cmask) 

19: append r ^ ( ■ Hl^ii-l)'^"'""''''-''^'''''"'''"'^ |Sforcc,l + Schoicc) (sforcc,2 + Schoicol 

20: t <^ NextFermi(i) 
21: end while 
22: return r 



2.4- Bosons 

As a short outlook, we want to illustrate how the developed methods can easily be 
adapted to bosonic systems as well. In quantum mechanics, bosons are subatomic parti- 
cles which obey Bose-Einstein statistics, like, for example, photons. For our purposes, we 
replace fermionic "orbitals" by bosonic "modes", which can be multiply occupied (i.e., 
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the Pauli exclusion principle no longer holds for bosons). That is, the bosonic analogue 
of a fermionic Slater determinant differs only by the unrestricted number of particles in 
each mode. The central observation of this subsection states that a bit-encoding (equiv- 
alent to Slater determinants) works for bosons as well. The idea is detailed in figure 5, 
where 0-bits serve as delimiters between modes. Lexicographical enumeration of bosonic 
states with a fixed total particle number N and number of modes m is accomplished via 
enumeration of the bit-encoded Slater determinants with (to -|- — 1) orbitals and N 
particles! That is, Algorithm 1 may be employed without modifications. 



1: 1 



2:2 



3: 



4:3 



I 

1 

I 



110 111 



I 



I 



Figure 5: Bit-encoding of a bosonic state. Blue numbers label modes, and 0-bits serve as delimiters 
between modes. The shown state consists of one boson in the 1st mode, two in the 2nd mode, zero in 
the 3rd and three in the 4th. 



3. Application to Transition Metal Atoms 

The application example is based on the series [15, 16, 17], in which [17] makes use 
of the FermiFab toolbox to calculate ground state approximations for transition metal 
atoms (employing so-called configuration-interaction (CI) methods). The underlying 
quantum mechanical (non-relativistic, Born-Oppenheimer) Hamiltonian H = Hq + Vee 
with 



i<j<i<jv 



governs atoms/ions with N electrons and nuclear charge Z. The two terms in Hq are the 
single-particle kinetic energy and nuclear potential, respectively, whereas the Coulomb 
operator Vee describes the pairwise inter-electron Coulomb repulsion. The Hamiltonian 
leaves the simultaneous eigenspaces of the well-known angular momentum, spin and 
parity ('LS') operators invariant, so calculating these eigenspaces first leads to a huge 
dimension reduction. Specifically, the FermiFab toolbox automates the LS-eigenspace 
computation by combining configurations (4) with Clebsch-Gordan coefiicients. We skip 
further details here; instead, for the purpose of this section, we provide two orthonormal 
LS-eigenstates of neutral Chromium {N = Z ^ 24) with symmetry level ^D: 

ipi -.— —j= ( \ 3do 3dm. Sdx 4s 4s 4dx ) — 1 3do 3dm Sdy 4s 4s Ady ) — 1 3do 3dz 3dx 4s 4s 4dy ) 



( 3d(, 3dz 3dy 4s 4s 4dj; ) + \/3 ( Sd^ 3dm Sd^ 4s 4s 4dy ) — -\/3 ( 3dz 3dm 3dy 4s 4s id^ 
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and 

V'2 := (-\/3/2|3do4s4s4p^4p^4dy) - V3/2|3do4s4s4p^4py4d^) 



. 21 

+ 2 \3dm4sla4p^ '^Py'^d^) + ^ j3d„4s4s4p2 4p:,4rfj;) - ^ |3d„ AsIsAp^ 4py Ad^) 

+ ^ \3d^ 4sls 4p^ Apy 4dy) + \3d^ 4s4s4p^4p^ 4d^) + Vs\3d:c4s4s4p^4py 4do) 

— \3dx 4s4s4p3 4py 4dm) — ^ \3dy 4s4s4p^ 4py 4dx) + V3\3dy 4s4s4p^ 4p^ 4do) 
+ \3dy4s4s4pz 4px 4dm) + \3dy4s4s 4p^ 4py 46?^) — 2|3d24s4s 4px 4py 4dm) 

+ ^ \ 'id^4s4s4p^4px 4d^) + ^ \ Mz4s4s4p^4py4dy) 

In this notation, ~ means spin down ^, otherwise up fi and the s,p,d subshell orbitals 
are labeled s, PzPxPy and do dz dm dx dy, respectively. The numbers 3 and 4 denote the 
third and fourth shell. Since all spin-orbitals up to 3p are fully occupied, they are not 
shown here for conciseness of notation. 

The following paragraph demonstrates how to translate the expectation value ('02 | Vee V'l) 
into a list of Coulomb integral symbols 

{ab\cd):= [ a*{xi)b{xi)- r c* {x2)d{x2) dxi dx2, (9) 

JrB \Xi-X2\ 

where a,b,c,d € L^(]R'^) are spatial orbitals and * denotes complex conjugation. As 
shown in (8), the essential step is the calculation of the 2-body reduced density matrix 
^\i>i) {■4>2\ ■ Using the Mathematica version of FermiFab, this is accomplished by the first 
line of the following code sample (see mathematica/RDMdemo .nb); the subsequent code 
just displays the result: 

p2 = RDM[l/ri, ll/2, 2] ; 

SlaterToPureString [orbnaines_ , s_] : = 
StringJoin @@ Flatten [ {B, " " } & /@ orbnames [[FermiToCoords [s] J ] [l ; ; -2J 

(* display RDM *) 

P2 /. IXI[ket_, bra_] " | "<> SlaterToPureString [allShells , ket] <> 
"><"<> SlaterToPureString [allShells , bra] <>"|" 

3d0 3dx><4pz 4pY | | 3d0 3dy><4pz 4px | | 3dz 3dT[i><4px 4py 



V210 V210 ^fW 



The FermiToCoords command converts any bit-encoded Slater determinant to a vector 
of integers enumerating the occupied orbitals. 

Since the Coulomb operator is independent of spin, we may effectively "trace out" 
the spin coordinate from the employed spin-orbitals. Specifically, consider single-particle 
wavefunctions 

Xt{x,s) ^ ^p^{x)a^{s), cc G M^ s G {t, 1} , i = l,...,4 
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which factor into the spatial part Lpi and spin part a^. Endowing particle i with coordi- 
nates {xi,Si), the antisymmetrized 2-body Slater determinants read 

\XiX3) = ^ {Vi{xi)ai{si) tfj{x2)aj{s2) - </Jj(a;i)aj(si) (pi(a;2)ai(s2)) ■ 
Plugged into the following equation for the Coulomb expectation value yields 

( Xi X2 I ] ^- rX3 X4 ) 

\ \xi-X2\ I 

= ((/3i(^3 I </'2'^4) {oi\ I as) {PL2 I Q!4) 

- ((/?i(y54 I ¥'29'3) ("1 I "4) ("2 I as) • 

Translating this equation to alternating spin up f and down J, orbitals (and taking sym- 
metries of (ah I cd) into account) is accomplished by the SpinTraceCoulomb command in 
the first line of the following code sample: 

pSj = SpinTraceCoulomb [ (List @@ P2) /• c_IXI[ket_, bra_] {c, {ket, bra}}]; 

BosonToString [orbnames_ , s_] := 
StringJoin @@ Flatten [ {tt, " " } & /@ orbnames IBosonToCoords [s] J ] [l ; ; -2]] 

(* display Coulomb integrals *) 

pSj /• Coulint [a_, b_] :-» "(" <> BosonToString [allShellsSpatial , a] <> 
" I " <> BosonToString [allShellsSpatial , b] <> ") " 

(3d0 4px|3dy 4pz) (3d0 4pY | 3dx 4pz) (3d0 4pz | 3dx 4pY) 

+ + 

V2IO V2IO V2IO 

(3d0 4pz|3dy 4px) (3dz 4px|3dm 4py) (3dz 4py | 3dm 4px) 

+ - 

V210 a/7o~ 

Note that spatial orbitals can appear twice within a Coulomb integral symbol, e.g., 
(aa\hc). Thus, a bosonic encoding of these spatial orbitals is used to accommodate 
multiple occurrences, and hence the BosonToCoords command. 

Concluding, we have obtained the desired list of Coulomb integral symbols, which 
may then be evaluated by inserting concrete functions into (9). 
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